#ifndef SYMMULKERNEL_CU
#define SYMMULKERNEL_CU
#include "symMul.h"
#include <cuda_runtime.h>

__global__ void symMulKernel(float* M, float* B)
{
    int blockRow = blockIdx.y;
    int blockCol = blockIdx.x;

    float* Bsub;
    Bsub = &B[N*BLOCK_SIZE*THREAD_SIZE*blockRow+BLOCK_SIZE*THREAD_SIZE*blockCol];

    float Bvalue[THREAD_SIZE][THREAD_SIZE];

    int row =threadIdx.y;
    int col =threadIdx.x;

    for (unsigned int i=0; i<THREAD_SIZE; i++){
        for (unsigned int j=0; j<THREAD_SIZE; j++){
            Bvalue[i][j] = 0;
        }
    }

    for (unsigned int k=0; k< GRID_SIZE; k++){
        float* Msub;
        Msub = &M[N*BLOCK_SIZE*THREAD_SIZE*blockRow+BLOCK_SIZE*THREAD_SIZE*k];
        float* Mtsub;
        Mtsub= &M[N*BLOCK_SIZE*THREAD_SIZE*blockCol+BLOCK_SIZE*THREAD_SIZE*k];

        __shared__ float Ms[BLOCK_SIZE*THREAD_SIZE][BLOCK_SIZE*THREAD_SIZE];
        __shared__ float Mts[BLOCK_SIZE*THREAD_SIZE][BLOCK_SIZE*THREAD_SIZE];

        for (unsigned int i=0; i<THREAD_SIZE; i++){
            for (unsigned int j=0; j<THREAD_SIZE; j++){
                unsigned int x=row*THREAD_SIZE+i;
                unsigned int y=col*THREAD_SIZE+j;
                Ms[x][y]= Msub[x*N+y];
                Mts[x][y]= Mtsub[x*N+y];

            }
        }

        __syncthreads();

        for (int e=0; e<BLOCK_SIZE*THREAD_SIZE; e++){
            for (unsigned int i=0; i<THREAD_SIZE; i++){
                for (unsigned int j=0; j<THREAD_SIZE; j++){
                    unsigned int x=row*THREAD_SIZE+i;
                    unsigned int y=col*THREAD_SIZE+j;
                    Bvalue[i][j] += Ms[x][e]*Mts[y][e];
                }
            }
        }

        __syncthreads();
    }

    for (unsigned int i=0; i<THREAD_SIZE; i++){
        for (unsigned int j=0; j<THREAD_SIZE; j++){
            unsigned int x=row*THREAD_SIZE+i;
            unsigned int y=col*THREAD_SIZE+j;
            Bsub[idx(x,y)] = Bvalue[i][j];
        }

    }
}

#endif
